{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "\n",
    "#My name: Joe Puccio\n",
    "#Collaborators: Fred Landis, Alan Wu, Zen Yang\n",
    "\n",
    "#Helper functions courtesy of Master McMillan\n",
    "def GlobalAlign(v, w, match, mismatch, gap):\n",
    "\ts = zeros((len(v)+1,len(w)+1), dtype=int32)\n",
    "\tb = zeros((len(v)+1,len(w)+1), dtype=int32)\n",
    "\n",
    "\tfor j in xrange(1,len(w)+1):\n",
    "\t\ts[0,j] = s[0,j-1] + gap\n",
    "\t\tb[0,j] = 1 \n",
    "\n",
    "\tfor i in xrange(1,len(v)+1):\n",
    "\t\ts[i,0] = s[i-1,0] + gap\n",
    "\t\tb[i,0] = 2\n",
    "\n",
    "\tfor i in xrange(1,len(v)+1):\n",
    "\t\tfor j in xrange(1,len(w)+1):\n",
    "\t\t\tdiag = s[i-1, j-1] + (match if (v[i-1] == w[j-1]) else mismatch)\n",
    "\t\t\tvskip = s[i-1, j] + gap\n",
    "\t\t\twskip = s[i, j-1] + gap\n",
    "\t\t\ts[i,j] = max(vskip, wskip, diag)\n",
    "\t\t\tif (s[i,j] == wskip):\n",
    "\t\t\t\tb[i,j] = 1\n",
    "\t\t\telif (s[i,j] == vskip):\n",
    "\t\t\t\tb[i,j] = 2\n",
    "\t\t\telse:\n",
    "\t\t\t\tb[i,j] = 3\n",
    "\n",
    "\treturn s,b \n",
    "\t#returns the score array and the backtrack array\n",
    "\n",
    "def GlobalAlignForSubstring(v, w, match, mismatch, gap):\n",
    "\ts = zeros((len(v)+1,len(w)+1), dtype=int32)\n",
    "\tb = zeros((len(v)+1,len(w)+1), dtype=int32)\n",
    "\n",
    "\t#for j in xrange(1,len(w)+1):\n",
    "\t#\ts[0,j] = s[0,j-1] + gap\n",
    "\t#\tb[0,j] = 1 \n",
    "\n",
    "\tfor i in xrange(1,len(v)+1):\n",
    "\t\ts[i,0] = s[i-1,0] + gap\n",
    "\t\tb[i,0] = 2\n",
    "\n",
    "\tfor i in xrange(1,len(v)+1):\n",
    "\t\tfor j in xrange(1,len(w)+1):\n",
    "\t\t\tdiag = s[i-1, j-1] + (match if (v[i-1] == w[j-1]) else mismatch)\n",
    "\t\t\tvskip = s[i-1, j] + gap\n",
    "\t\t\twskip = s[i, j-1] + gap\n",
    "\t\t\ts[i,j] = max(vskip, wskip, diag)\n",
    "\t\t\tif (s[i,j] == wskip):\n",
    "\t\t\t\tb[i,j] = 1\n",
    "\t\t\telif (s[i,j] == vskip):\n",
    "\t\t\t\tb[i,j] = 2\n",
    "\t\t\telse:\n",
    "\t\t\t\tb[i,j] = 3\n",
    "\n",
    "\treturn s,b \n",
    "\t#returns the score array and the backtrack array\n",
    "\n",
    "def alignWNonRecursive(b,w,i,j,local=False):\n",
    "\twhile (b[i,j] != 0) and local:\n",
    "\t\tif (b[i,j] == 3):\n",
    "\t\t\ti = i-1\n",
    "\t\t\tj = j-1\n",
    "\t\telif (b[i,j] == 2):\n",
    "\t\t\ti = i-1\n",
    "\t\telif (b[i,j] == 1):\n",
    "\t\t\tj = j-1\n",
    "\treturn j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Problem 1\n",
    "\n",
    "#It seems like this is just asking us to do Local Alignment between the fish and human sequence, and then the \n",
    "#difference from normal local alignment is that we're going to choose the alignment that starts in the first \n",
    "#colmun and ends in the last column (assuming human gene spans top row), because we need to make sure that the \n",
    "#resulting substring from the fish sequence actually has the start and end matching of the human sequence. \n",
    "\n",
    "with open(\"humanGene.seq\", 'r') as sequenceFile:\n",
    "\thGene = sequenceFile.read()\n",
    "with open(\"mouseGene.seq\", 'r') as sequenceFile:\n",
    "\tmGene = sequenceFile.read()\n",
    "with open(\"yeastGene.seq\", 'r') as sequenceFile:\n",
    "\tyGene = sequenceFile.read()\n",
    "with open(\"zfishGene.seq\", 'r') as sequenceFile:\n",
    "\tfGene = sequenceFile.read()\n",
    "\n",
    "scores, backtrack = GlobalAlignForSubstring(hGene,fGene,1,0,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[    0     0     0 ...,     0     0     0]\n",
      " [   -1     0     0 ...,     1     1     1]\n",
      " [   -2    -1     1 ...,     1     1     1]\n",
      " ..., \n",
      " [-2725 -2723 -2721 ...,  1609  1608  1607]\n",
      " [-2726 -2724 -2722 ...,  1611  1610  1609]\n",
      " [-2727 -2725 -2723 ...,  1613  1612  1611]]\n"
     ]
    }
   ],
   "source": [
    "print scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0 0 ..., 0 0 0]\n",
      " [2 3 3 ..., 3 3 3]\n",
      " [2 2 3 ..., 3 3 3]\n",
      " ..., \n",
      " [2 2 2 ..., 1 1 1]\n",
      " [2 2 2 ..., 1 1 1]\n",
      " [2 2 2 ..., 1 1 1]]\n"
     ]
    }
   ],
   "source": [
    "print backtrack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3038"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "argmax(scores[len(hGene)]) \n",
    "#this gives us the end index we want for fish gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "317\n"
     ]
    }
   ],
   "source": [
    "print alignWNonRecursive(backtrack,fGene,len(hGene),3038,True) \n",
    "#this gives us the start index we want for the fish gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#So our start index of the fish gene is 317, and our end index is 3038. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_all_substrings_return_dictionary_with_initial_count(string):\n",
    "  length = len(string)\n",
    "  adict = {}\n",
    "  for i in xrange(length):\n",
    "    for j in xrange(i,length):\n",
    "      adict[string[i:j + 1]]=1 \n",
    "  return adict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGAAGAGGA\n"
     ]
    }
   ],
   "source": [
    "substringsOfYeast = get_all_substrings_return_dictionary_with_initial_count(yGene)\n",
    "substringsOfHuman = get_all_substrings_return_dictionary_with_initial_count(hGene)\n",
    "substringsOfMouse = get_all_substrings_return_dictionary_with_initial_count(mGene)\n",
    "substringsOfFish = get_all_substrings_return_dictionary_with_initial_count(fGene)\n",
    "\n",
    "#we save runtime if we only look at values that have been in all of them so far\n",
    "\n",
    "substringsInYH = []\n",
    "for substring in substringsOfYeast:\n",
    "\tif substring in substringsOfHuman: \n",
    "\t\tsubstringsInYH.append(substring)\n",
    "#print substringsInYH\n",
    "\n",
    "substringsInYHM = []\n",
    "for substring in substringsInYH:\n",
    "\tif substring in substringsOfMouse: \n",
    "\t\tsubstringsInYHM.append(substring)\n",
    "#print substringsInYHM\n",
    "\n",
    "substringsInAll= []\n",
    "for substring in substringsInYHM:\n",
    "\tif substring in substringsOfFish: \n",
    "\t\tsubstringsInAll.append(substring)\n",
    "\n",
    "print max(substringsInAll, key=len) \n",
    "#this prints the longest common substring in all of the gene sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Yeast and Human:  320\n",
      "Yeast and Mouse:  321\n",
      "Yeast and Fish:  306\n",
      "Fish and Human:  1723\n",
      "Fish and Mouse:  1722\n",
      "Mouse and Human:  2203\n"
     ]
    }
   ],
   "source": [
    "#Problem 3\n",
    "\n",
    "#Take the output of problem 1 in, because that's our last normalized sequence. \n",
    "fGene = fGene[317:3038]\n",
    "\n",
    "score, back = GlobalAlign(yGene, hGene,1,0,-1)\n",
    "print \"Yeast and Human: \",score[len(yGene)][len(hGene)]\n",
    "\n",
    "score, back = GlobalAlign(yGene, mGene,1,0,-1)\n",
    "print \"Yeast and Mouse: \",score[len(yGene)][len(mGene)]\n",
    "\n",
    "score, back = GlobalAlign(yGene, fGene,1,0,-1)\n",
    "print \"Yeast and Fish: \",score[len(yGene)][len(fGene)]\n",
    "\n",
    "score, back = GlobalAlign(fGene, hGene,1,0,-1)\n",
    "print \"Fish and Human: \",score[len(fGene)][len(hGene)]\n",
    "\n",
    "score, back = GlobalAlign(fGene, mGene,1,0,-1)\n",
    "print \"Fish and Mouse: \",score[len(fGene)][len(mGene)]\n",
    "\n",
    "score, back = GlobalAlign(mGene, hGene,1,0,-1)\n",
    "print \"Mouse and Human: \",score[len(mGene)][len(hGene)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}